################################################
#	Spry and Nunnally (2025)
# Replication Data for "Modern American Dilemma"
# Linear Models
################################################

# Front Matter ------------------------------------------------------------

rm(list=ls())		# clear objects in memory

setwd("/Users/amberspry/Dropbox/Research/Data/2016 CMPS/Analysis")
# Data imported from 2016 CMPS

source("2026_sprynunnally_source.R")

# Standard Error Function -------------------------------------------------

se_mean <- function(x){
  x_star <- x
  x_star <- x_star[!is.na(x)]
  se <- sd(x_star)/sqrt(length(x_star))
}

# Diagnostic Plots --------------------------------------------------------

# Descriptive Stats

# Race Tables
wtd.table(sub.bw$dummy.white, weights = sub.bw$weight)

# American ID by race tables
tab.amer <- questionr::wtd.table(sub.bw$dummy.black, sub.bw$amer, weights = sub.bw$weight)
prop.table(tab.amer,1)
df.amer <- data.frame(NotAtAll = c(4, 5), NotVery = c(6, 10), Somewhat = c(24, 27), Very = c(66, 58), Race = c("White", "Black"))
m.amer <- melt(df.amer, id.vars = 'Race')

ggplot(m.amer, aes(fill=variable, y=value, x=Race)) +
  geom_bar(position = "dodge", stat = "identity") +
  labs(x="Race", y="Percentage of Respondents", title = "Importance of American Identity by Racial Group") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_grey(name="Importance of American ID",
                  breaks=c("NotAtAll", "NotVery", "Somewhat", "Very"),
                  labels=c("Not At All", "Not Very", "Somewhat", "Very"))


# BLM Support by race tables

tab.blm <- questionr::wtd.table(sub.bw$dummy.black, sub.bw$blm, weights = sub.bw$weight)
prop.table(tab.blm,1)
df.blm <- data.frame(One = c(31, 5), Two = c(13, 3), Three = c(30, 22), Four = c(15, 29), Five = c(11, 41), Race = c("White", "Black"))
m.blm <- melt(df.blm, id.vars = 'Race')

ggplot(m.blm, aes(fill=variable, y=value, x=Race)) +
  geom_bar(position = "dodge", stat = "identity") +
  labs(x="Race", y="Percentage of Respondents", title = "Attitudes Toward Black Lives Matter by Racial Group") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_grey(name="Attitudes\nToward Black Lives Matter",
                  breaks=c("One", "Two", "Three", "Four", "Five"),
                  labels=c("Strongly Oppose", "Somewhat Oppose", "Neither Support nor Oppose", "Somewhat Support", "Strongly Support"))

# Linear Models for Manuscript --------------------------------------------

# Conditional OLS on white and black subsets within same dataset
# THESE ARE THE TABLES TO USE

# Model 1: Multivariate (Democrat + Independents) --------------------------------
# Multivariate analysis (Democrat and Independents Separate Variables)
# USE THIS MODEL

white.blm.sc.lmb <- lm(blm ~ amer
                       + female
                       + age
                       + dparty
                       + iparty
                       + protest
                       + bhood
                       + income
                       , data = subset(sub.bw, s2_1==1), weights = weight)
summary(white.blm.sc.lmb)

black.blm.sc.lmb <- lm(blm ~ amer
                       + female
                       + age
                       + dparty
                       + iparty
                       + protest
                       + bhood
                       + income
                       , data = subset(sub.bw, s2_3==1), weights = weight)
summary(black.blm.sc.lmb)

stargazer(white.blm.sc.lmb, black.blm.sc.lmb, title = "Results", align = TRUE,
          column.labels = c("White", "Black"),
          dep.var.labels = c("Attitudes Toward Black Lives Matter"),
          covariate.labels = c("American Identity", "Female", "Age", "Democrat", "Independent", "Support for Protest", "Neighborhood Percent Black", "Income"),
          omit.stat = c("ser", "f"),
          digits = 2,
          no.space = TRUE)

plot_summs(black.blm.sc.lmb, white.blm.sc.lmb, scale=TRUE,
           coefs = c("American"="amer", "Female"="female", "Age"="age", "Democrat"="dparty", "Independent"="iparty", "Support for Protest"="protest", "Neighborhood Percent Black"="bhood", "Income"="income"),
           model.names = c("Black", "White"))


# Model 2A: OLS, added controls, split by race --------------------
# OLS Regression plus additional controls (Democrat and Independents separate variables)

white.blm.sc.lm2a <- lm(blm ~ amer
                        + female
                        + age
                        + dparty
                        + iparty
                        + protest
                        + bhood
                        + income
                        + lf
                        + havesay
                        + help1
                        + efficacy
                        + cell
                        , data = subset(sub.bw, s2_1==1), weights = weight)
summary(white.blm.sc.lm2a)

black.blm.sc.lm2a <- lm(blm ~ amer
                        + female
                        + age
                        + dparty
                        + iparty
                        + protest
                        + bhood
                        + income
                        + lf
                        + havesay
                        + help1
                        + efficacy
                        + cell
                        , data = subset(sub.bw, s2_3==1), weights = weight)

summary(black.blm.sc.lm2a)

stargazer(white.blm.sc.lm2a, black.blm.sc.lm2a, title = "Results", align = TRUE,
          column.labels = c("White", "Black"),
          dep.var.labels = c("Attitudes Toward Black Lives Matter"),
          covariate.labels = c("American Identity", "Female", "Age", "Democrat", "Independent", "Support for Protest", "Neighborhood Percent Black", "Household Income", "Linked Fate", "Racial Group Members Have a Say", "Public Officials Help My Group", "What People Like Me Think", "Recorded Police Misconduct"),
          omit.stat = c("ser", "f"),
          digits = 2,
          no.space = TRUE)
library(jtools)
plot_summs(black.blm.sc.lm2a, white.blm.sc.lm2a, scale=TRUE,
           coefs = c("American"="amer","Female"="female", "Age"="age", "Democrat"="dparty", "Independent"="iparty", "Support for Protest"="protest", "Neighborhood Percent Black"="bhood", "Income"="income", "Linked Fate"="lf", "Racial Group Members Have a Say"="havesay", "Public Officials Help My Group"="help1", "What People Like Me Think"="efficacy", "Recorded Police Misconduct"="cell"),
           model.names = c("Black", "White"),
           colors = c("black", "gray50"))

# Model 2B: OLS, added controls, American x Race (Dem + Indep) -----------------------
# OLS Regression plus additional controls (Democrat and Independents separate variables)
# Same as Model 2 but with interaction instead of split by race

int.blm.lm <- lm(blm ~ amer*dummy.white
                 + female
                 + age
                 + dparty
                 + iparty
                 + protest
                 + bhood
                 + income
                 + lf
                 + havesay
                 + help1
                 + efficacy
                 + cell
                 , data = subset(sub.bw), weights = weight)
summary(int.blm.lm)

stargazer(int.blm.lm, title = "Results - Race * American Interaction Model", align = TRUE,
          dep.var.labels = c("Attitudes Toward Black Lives Matter"),
          covariate.labels = c("American Identity", "White", "Female", "Age", "Democrat", "Independent", "Support for Protest", "Neighborhood Percent Black", "Household Income", "Linked Fate", "People Like Me Have a Say", "Public Officials Help My Group", "What People Like Me Think", "Recorded Police Misconduct", "American * White"),
          omit.stat = c("ser", "f"),
          digits = 2,
          no.space = TRUE)

# Model 3: Amer x LF, Amer x Partisanship, Amer x Gender Interactions (R Holdout) -------
# LF, Partisanship, Gender Interactions

white.blm.sc.lf3 <- lm(blm ~ amer*lf
                       + amer*dparty
                       + amer*iparty
                       + amer*female
                       + age
                       + protest
                       + bhood
                       + income
                       + havesay
                       + help1
                       + efficacy
                       + cell
                       , data = subset(sub.bw, s2_1==1), weights = weight)
summary(white.blm.sc.lf3)

black.blm.sc.lf3 <- lm(blm ~ amer*lf
                       + amer*dparty
                       + amer*iparty
                       + amer*female
                       + age
                       + protest
                       + bhood
                       + income
                       + havesay
                       + help1
                       + efficacy
                       + cell
                       , data = subset(sub.bw, s2_3==1), weights = weight)
summary(black.blm.sc.lf3)

stargazer(white.blm.sc.lf3, black.blm.sc.lf3, title = "Results", align = TRUE,
          column.labels = c("White", "Black"),
          dep.var.labels = c("Attitudes Toward Black Lives Matter"),
          covariate.labels = c("American Identity", "Linked Fate", "Democrat", "Independent", "Female", "Age", "Support for Protest", "Neighborhood Percent Black", "Household Income", "Racial Group Members Have a Say", "Public Officials Help My Group", "What People Like Me Think", "Recorded Police Misconduct", "American * Linked Fate", "American * Democrat", "American * Independent", "American * Female"),
          omit.stat = c("ser", "f"),
          digits = 2,
          no.space = TRUE)

plot_summs(black.blm.sc.lf3, white.blm.sc.lf3, scale=TRUE,
           coefs = c("American*Linked Fate"="amer:lf","American*Democrat","American*Independent", "American*Female"="female", "American"="amer","Linked fate"="lf", "Female"="female", "Age"="age", "Democrat"="dparty", "Independent"="iparty", "Support for Protest"="protest", "Neighborhood Percent Black"="bhood", "Income"="income", "Racial Group Members Have a Say"="havesay", "Public Officials Help My Group"="help1", "What People Like Me Think"="efficacy", "Recorded Police Misconduct"="cell"),
           model.names = c("Black", "White"))


# Model 4: Amer x LF, Amer x Partisanship, Amer x Gender Interaction (D Holdout) --------
# LF, Partisanship, Gender Interactions

white.blm.sc.lf4 <- lm(blm ~ amer*lf
                       + amer*rparty
                       + amer*iparty
                       + amer*female
                       + age
                       + protest
                       + bhood
                       + income
                       + havesay
                       + help1
                       + efficacy
                       + cell
                       , data = subset(sub.bw, s2_1==1), weights = weight)
summary(white.blm.sc.lf4)

black.blm.sc.lf4 <- lm(blm ~ amer*lf
                       + amer*rparty
                       + amer*iparty
                       + amer*female
                       + age
                       + protest
                       + bhood
                       + income
                       + havesay
                       + help1
                       + efficacy
                       + cell
                       , data = subset(sub.bw, s2_3==1), weights = weight)
summary(black.blm.sc.lf4)

stargazer(white.blm.sc.lf4, black.blm.sc.lf4, title = "Results", align = TRUE,
          column.labels = c("White", "Black"),
          dep.var.labels = c("Attitudes Toward Black Lives Matter"),
          covariate.labels = c("American Identity", "Linked Fate", "Republican", "Independent", "Female", "Age", "Support for Protest", "Neighborhood Percent Black", "Household Income", "Racial Group Members Have a Say", "Public Officials Help My Group", "What People Like Me Think", "Recorded Police Misconduct", "American * Linked Fate", "American * Republican", "American * Independent", "American * Female"),
          omit.stat = c("ser", "f"),
          digits = 2,
          no.space = TRUE)


# Plot the Interaction Effects --------------------------------------------
# Using Model 2 (American * Race interaction)

library(ggeffects)

pred_values <- ggpredict(int.blm.lm, terms = c("amer", "dummy.white"))

int.plot1 <- plot(pred_values) +
  labs(title = "Predicted Values of Support for Black Lives Matter: Model 2",
       x = "American Identity",
       y = "Support for Black Lives Matter") +
  scale_color_manual(values = c("black", "grey60"),
                     labels = c("Black", "White"),
                     name = "Race") +
  scale_fill_manual(values = c("black", "grey60"),
                    labels = c("Black", "White"),
                    name = "Race") +
  scale_y_continuous(limits = c(1,5)) +
  theme_bw()

print(int.plot1)
